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Relativistic dissipative fluid dynamics is a common tool to describe the space-time evolution of the 
strongly interacting matter created in ultrarelativistic heavy-ion collisions. For a proper comparison 
to experimental data, fluid-dynamical calculations have to be performed on an event-by-event basis. 
Therefore, fluid dynamics should be able to reproduce, not only the event-averaged momentum 
anisotropics, (v n ), but also their distributions. In this paper, we investigate the event-by-event 
distributions of the initial-state and momentum anisotropics e n and v n , and their correlations. We 
demonstrate that the event-by-event distributions of relative v n fluctuations are almost equal to the 
event-by-event distributions of corresponding e n fluctuations, allowing experimental determination 
of the relative anisotropy fluctuations of the initial state. Furthermore, the correlation c(v2,Va) turns 
out to be sensitive to the viscosity of the fluid providing an additional constraint to the properties 
of the strongly interacting matter. 



I. INTRODUCTION 



Relativistic dissipative fluid dynamics is the most widely employed model to describe the space-time evolution of 
the quark-gluon plasma (QGP) created in ultrarelativistic heavy-ion collisions. It was the success of fluid-dynamical 
models in describing the large azimuthal momentum anisotropics observed in heavy-ion collisions that led to our 
current picture of the QGP, as a strongly interacting fluid with one of the smallest shear viscosity to entropy density 
ratios, rj/s, ever observed 

The azimuthal momentum anisotropy is characterized in terms of the coefficients v n of the Fourier expansion of the 
single particle azimuthal distribution: 
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where ip n is the event- plane angle, ip n = (1/n) arctan((pTsinn0) / (pr cos n<p)), and <p is the transverse momentum 
azimuthal angle. One of the main features of the fluid-dynamical description of the expansion is that the anisotropy 
originates from the azimuthal anisotropy of the initial density profile. In the literature this initial anisotropy is 
quantified in terms of coefficients e TOj „: 
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where e is the energy density, r 2 = x 2 + y 2 , <p is now the spatial azimuthal angle, and ^ m .n is the participant angle, 
defined as 
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In the following we will concentrate on the anisotropics £2,ni an d use a shorthand e„ = £2,n- 

In fluid-dynamical calculations, a linear relation between vi and £2 was found, i.e., 1*2 oc e 2 0- The proportionality 
coefficient was shown to depend, not only on the properties of the fluid such as the equation of state and viscosity, 
but also on the initial density, freeze-out temperature, and resonance content of the late hadronic state 043l- 

The initial conditions in these fluid dynamical calculations were always smooth, constructed as an average over 
infinitely many individual collisions of the particular centrality. It was thought that the use of this kind of averaged 
initial conditions would lead to a good description of observables which were averaged over many events. In other 
words, the anisotropy v 2 computed using the event-averaged initial condition was expected to be equal to the (v 2 
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observed in the collisions, where (. . .) cv corresponds to an average over all the events of the corresponding centrality 
class. 

Recently, it has been realized that, in order to obtain a proper comparison with experimental data, fluid-dynamical 
calculations have to be performed on an event- by- event basis. This was first pointed out by Kodama et al. [5j almost 
ten years ago. However, this view only become widely accepted years later, after the work of Alver and Roland [f|. 
They showed that 63 and, consequently, 1*3 are non-zero in a single event. This is in contrast to the traditionally 
used event-averaged initial conditions in fluid-dynamical models which had zero 63 and V3. Furthermore, Alver and 
Roland demonstrated that such finite value of U3 can be observed in heavy-ion collisions. This finding made (w3) ov as 
important observable as (v2) cv for probing the properties of the dense matter formed in heavy-ion collisions, and led 
to several works studying the behavior of observables in an event-by-event fluid-dynamical description 0, H| . 

On the other hand, if fluid dynamics can be applied to describe individual ultrarelativistic heavy-ion collisions, it 
must be able to describe v n in every collision, not only the average (v n ) cv . Therefore it must be able to reproduce the 
distribution V(v n ) of v„ in an ensemble of events too. To confirm the applicability of fluid dynamics to describe the 
expansion stage of heavy-ion collisions, it is thus not enough to check whether the event-averaged values of v n agree 
with the data, but one must also check whether their distributions, V (v n ), match what is experimentally observed. 
Recently, the distributions of V2, i>3, and V4 were measured at the LHC by the ATLAS collaboration Q. Also, the 
first fluid-dynamical calculations of these distributions were performed by Gale et al. [Io| . 

In this paper, we study the event-by-event probability distribution of the Fourier coefficients v n , P(v n ), and how 
they are correlated with the initial state anisotropics e n event-by-event. The goal of this paper is not to attempt 
a comparison with experimental data, but to explore how these distributions and correlations are affected by the 
fluid viscosity and initialization of the system. In this way, it will be possible to understand what can be learned by 
measuring such event-by-event distributions. 

In the following we explain our fluid dynamical model in section UH and show our results in section Mil Subsection 
IIII Al is dedicated to an analysis of the event- by- event correlation between initial condition and flow anisotropy, while 
in subsections IIII Bl and IlLICl we show our results for probability distributions of scaled anisotropy 6v n , P(6v n ), 
and linear correlation coefficients c(v n ,v m ), respectively. In section ITVl we summarize our findings and make our 
conclusions. 



II. MODEL 

To generate the initial states event-by-event, we use a Monte-Carlo Glauber model as implemented in Ref. Q. In 
this model, nucleons are distributed into nuclei according to Woods-Saxon distribution. NN-correlations and finite 
size effects are neglected since they have a negligible effect on the anisotropy coefficients [Tl[. In an event with a 
given impact parameter, nucleons from different nuclei are assumed to collide when their transverse distance d is small 
enough, i.e., when d 2 < <jnn/it. 

We consider two initial conditions, in which the initial entropy density, s, at tq — 1 fm, is evaluated as 



part ,bi] 



s(x,y) = W exp{-[(x- Xi ) 2 + (y- yi ) 2 ]/(2a 2 )}, (4) 

i=l 

where Xi and yi are the spatial coordinates of either wounded nucleons (initial condition sWN) or binary collisions 
(initial condition sBC), given by the Monte-Carlo Glauber model. W is a normalization constant fixed to provide the 
observed multiplicity and a = 0.8 fm is the spatial scale of a wounded nucleon or a binary collision. The centrality 
classes are determined according to the number of binary collisions (for initial condition sBC) or the number of 
participants (for initial condition sWN). The initial fluid velocity and shear-stress tensor are set to zero and we 
neglect the effects of bulk viscosity. 



For the fluid-dynamical evolution, we use the model previously employed in Ref. 12[ . We describe the time evolution 
of the fluid in the central rapidity region assuming boost invariancc and a zero baryochemical potential. The equations 
of motion are given by the conservation laws for energy and momentum: 

d„T^ = 0, (5) 

where T^ v = (e + p) u^u v — gVp + ir 1 - 1 " , with e, p, u^, and tt^ being the energy density, the thermodynamic pressure, 
the fluid 4-velocity, and the shear-stress tensor, respectively. We use the lattice QCD and hadron resonance gas 
based equation of state s95p-PCE-vl [l| with chemical freeze-out at temperature T c /j em = 150 MeV. The evolution 



equation of the shear-stress tensor is given by transient relativistic fluid dynamics [I4j, |15| : 
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where rj is the shear viscosity coefficient, D = u^d^ is the co-moving time derivative, = A'^'id a u 13 is the shear 
tensor, 6 = is the expansion rate, and A^ = (A£A;g + A^AJ, 1 - 2/SA^A^ /2, with A^ = g» v -u^u v . The 

transport coefficients of the non-linear terms on the right-hand side of the Eq. [S] were taken in the massless limit, 
in the f 4- moment approximation, and the relaxation time was assumed to be r w = 5?y/ (e + P) [TEl . Ilq ]. Here, we 
have not included the nonlinear terms related to the vorticity tensor. Note that the last two terms in Eq. ([6]) were 
not included in our previous studies [j^ . While such terms can have a significant effect on many observables, they 
are not relevant for the results discussed in this paper. We shall leave a detailed investigation of the effect of such 
terms to a future work. The equations of motion were solved numerically using the SHASTA algorithm, whereas the 
evolution equations for shear stress (Eq. [5]) were solved using simple finite differencing scheme. For more details see 
Refs. [H,[i|. 

The hadron spectra are calculated with the Cooper- Frye freeze-out procedure [l8j using the decoupling temperature 
Tf = 1 00 McV, which was shown to give reasonable agreement with both the py-spectrum and (1)2) ev for pions at 
RHIC when a temperature-dependent rj/s was used, see Refs. [T^ . In this work, we use constant values of viscosity, 
rj/s = and 0.16. Nevertheless, the pr-spectrum and (^2) ev remain close to what is actually observed at RHIC. Since 
our main purpose is not the comparison to experimental observables, we adjusted only the initial entropy density 
to fit the observed multiplicity, but kept all the other parameters unchanged. Finally, we use Israel and Stewart's 
I4-moment ansatz for the dissipative correction to the local equilibrium distribution function, 

M v 

where foi = {exp [(u^pf — jii) /T] ±1} 1 is the local equilibrium distribution function, with the index i indicating 
different hadron species and the four-momentum of the corresponding hadron. After calculating the thermal 
spectra, we include the contribution from all 2- and 3- particle decays of unstable resonances up to I.I GeV mass. 

It should be noted that because we do not generate particle ensembles at any point we always know the direction 
of the event plane and the magnitude of v n exactly. Experimentally, one measures a finite number of particles, which 
smears the observed distribution of v n . However, the final experimental result for the v n distributions undergoes an 
unfolding procedure that is supposed to remove such a smearing Q. Therefore, for a comparison with data, one can 
use the particle distributions computed with fluid dynamics without generating an ensemble of particles. A more 
detailed way would be to generate the particle ensembles and apply the same complicated unfolding procedure used 
by the experimentalist to obtain the v n distribution, but this procedure would be an unnecessary complication for 
the purpose of this work. 



III. RESULTS 



In this work we consider Au+Au collisions at ^snn = 200 A GeV. All the results shown in this paper are for 
positively charged pions. For each centrality class a total of 2000 events were computed. The Fourier coefficients and 
the initial-state anisotropics were calculated according to Eqs. (H), an d ©, respectively. In the following, we consider 
two constant values for the shear viscosity to entropy density ratio, rj/s = and 0.16. 



A. Correlations 



As mentioned in the Introduction, it has been known for a long time that the event averaged v 2 , and the eccentricity 
of the averaged initial state, £2 are a ppr oximately linearly related 0. Similar relation has been found for £3 and the 
average ^3 but not for £4 and V4 (l9l. |20| . Here we study whether similar relations hold event-by-event by evaluating 
the linear correlation between the harmonics v n and e n . We use the linear correlation coefficient, 

c(a, b ) = ( {a - {a) - nb - {bU ) , (8) 

\ °aOb / ev 

where <r a is the standard deviation of the quantity a. This correlation function is 1 (—1) if a and b arc linearly 
(anti-linearly) correlated and zero in the absence of linear correlation. 

The 2-dimensional histograms in Figs. [THU show the correlations between £2 and V2, £3 and V3, and £4 and W4, 
respectively, for the 20 — 30 % centrality class. To study the effect of both viscosity and initialization on these 
correlations, we show the correlations in three different cases: (a) sBC initialization with rj/s = 0, (b) sBC initialization 
with rj/s = 0.16, and (c) sWN initialization with rj/s = 0.16. 
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FIG. 1. 62 and V2 of pions in the 20 — 30 % centrality class using different initializations and viscosities, a) sBC and 77/s — 0, 
b) sBC and 77/s = 0.16 and c) sWN and 77/s = 0.16. 






FIG. 3. e 4 and 1)4 of pions in the 20 — 30 % centrality class using different initializations and viscosities, a) sBC and t]/s — 0, 
b) sBC and 77/s = 0.16 and c) sWN and 77/s = 0.16. 



As can be seen in these figures, the V2 and V3 coefficients display a strong linear correlation to their corresponding 
initial-state coefficients for all cases considered. This is confirmed by the values of the linear correlation coefficient 
c (v2, £2) ~ c (vs, 63) ~ 1, as shown in the Figures (top left corner). As for any two variables we can write 

v n = C n e n + 6 n , (9) 

where C n = (v n ) QV / (en) ev , and consequently, (<5„) cv = 0. The values of C n are shown in Figs. Q] and [2j For n = 2 
a linear relation, 1*2 = 6*262, is approximately satisfied event-by-event with only ~ 10% deviations from this relation 
at a given 62. On the other hand, an event-by-event linear relation between V3 and 63 is not satisfied well, with ^3 
deviating ~ 50% from tj 3 = 6363 at a given 63. 

In all the cases considered above, there is basically no linear correlation between 64 and V4, see Fig. [3] At least 
one reason for this behavior is that there is also correlation between e| and V4, which can be of the same order or 
larger than c(e 4 ,u 4 ): c(e|,u 4 ) = 0.40 (sBC, 77/s = ), c(ei,u 4 ) = 0.69 (sBC, 77/s = 0.16) and, c{e 2 2 ,v A ) = 0.46 (sWN, 
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rj/s = 0.16). This is a non-linear effect triggered by the coupling between two different Fourier coefficients. 



21|. 



and n = 4, and a linear combination of these two components was found to be a good estimator for 

As expected, the proportionality coefficients C n are sensitive to the value of the shear viscosity. This can be seen 
by comparing Figs. [T^, and [lb (n = 2), Figs. Ot and [2b (n = 3), and Figs. [3Jt and [5b (n = 4). In general, the higher 
Fourier coefficients are expected to be more sensitive to the viscosity [19j . This is also the case in our calculations, 
and is confirmed by comparing the relative changes in the coefficients C2, C3 and C4. 

Note that the proportionality constants C n do not depend only on the intrinsic properties of the fluid, but also on 
the initial conditions. Again something to be expected, since in the calculations done using averaged initial conditions, 
the precise value of the proportionality depended on many details as discussed in the Introduction. 
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FIG. 4. £2 and V2 in the — 5% centrality class using different initializations and viscosities, a) sBC and rj/s — 0, b) sBC and 
rj/s = 0.16 and c) sWN and rj/s = 0.16. 




FIG. 5. £3 and V3 in the — 5% centrality class using different initializations and viscosities, a) sBC and rj/s — 0, b) sBC and 
ij/s = 0.16 and c) sWN and rj/s = 0.16. 





FIG. 6. £4 and v± in the — 5% centrality class using different initializations and viscosities, a) sBC and rj/s — 0, b) sBC and 
rj/s = 0.16 and c) sWN and rj/s = 0.16. 



In Figs. S] and [5] we show the 2-dimcnsional histograms of €2 and i>2 and of £3 and V3, respectively, in the — 5 % 
centrality class. We plot the same cases considered above: (a) sBC initialization with rj/s = 0, (b) sBC initialization 
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with rj/s — 0.16, and (c) sWN initialization with rj/s = 0.16. For n — 2 and n = 3 the linear correlation is still 
valid. Also, the effect of shear viscosity and initialization on C n remain qualitatively the same. On the other hand, 
in Fig. |6] the correlation between £4 and V4 in central collisions is drastically different from the correlation in the 
20 — 30 % centrality class. In the — 5 % centrality class the linear correlation coefficient 0(64,^4) becomes much 
closer to 1 when compared to the peripheral case. It can be as large as ~ 0.81 obtained for the sWN initialization 
with rj/s = 0.16. This behavior is expected since in Ref. [2~i| it was shown that €4 becomes a better estimator for V4 
in central collisions. 



B. distributions of v. 



So far the event-averaged values of v n have been extensively studied. In order to observe what can be learned by 
looking at v„ probability distributions, it is convenient to remove the average from the distributions, and study the 
relative fluctuations using the scaled variables 

r V n — \V n )cv , r £n ~ \ e n/cv /, n \ 

5v n = — - — r , and 6e n = — - — . (10) 

In this way changes in the probability distributions due to changes in the average values are removed. 

It was shown in the previous subsection that v n and e n have a strong linear correlation for n = 2 and 3. As discussed 
in the Appendix, if two variables are linearly correlated, and (d) =0, the variances of the relative distributions are 
equal. Since viscosity has only a small effect on the correlations of v„ and e„, we expect that a% m erf , independent of 
viscosity. In such a case the information about the fluid response to the initial geometry is contained in the coefficients 
C n controlling the average (v n ) cv , while the relative fluctuations of v n originate from the relative fluctuations of e„ 
and do not depend on viscosity of the fluid. 




FIG. 7. Probability distributions: a) P(Sv 2 ) and P(Se 2 ), b) P(Sv 3 ) and P(6e 3 ), and c) P(5vi) and P(5e 4 ), in the 20 - 30 % 
centrality class with sBC initialization and two different values of rj/s, rj/s — and rj/s = 0.16. 




FIG. 8. Probability distributions: a) P(Sv 2 ) and P(Se 2 ), b) P(Sv 3 ) and P(Se 3 ), and c) P(8v 4 ) and P(5e 4 ), in the 20 - 30 % 
centrality class with rj/s = 0.16 and two different initial conditions, sBC and sWN. 

To test this assumption, and to sec whether similarity extends beyond variances, we plot the probability distributions 
P(8v n ) and P(Se n ) in 20 — 30 % centrality class in Fig. [7] using the sBC initialization and two values of viscosity, 
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rj/s = and 0.16. As seen in panel a) for n = 2 and in panel b) for n = 3, not only the variances are similar, but the 
entire distributions are almost identical. Fig. [7£ depicts the relative distributions for n = 4, and surprisingly they are 
very similar even if 1)4 and £4 are not linearly correlated. There are deviations only at the tail of the distribution. As 
discussed in the Appendix, lack of correlation leads to a large spread of possible values of 1*4 at given £4 which tends 
to make Sv^ distribution wider, but this effect can be canceled by other terms. How these terms arise is an interesting 
question beyond the scope of this paper. 

To check whether the similarity of P(Sv n ) and P(Se n ) is only a coincidence based on the sBC initialization, we show 
in Fig. [5] the distributions using both sBC and sWN initialization, but using only one value of viscosity, rj/s = 0.16. 
Again, panels a, b, and c depict cases with n = 2, 3, and 4, respectively. P(8v n ) and P(Se n ) are almost equal for both 
initializations. The distributions P(8e n ) are also similar for both initializations, but this is because both Glauber-type 
initializations give rise to the same relative fluctuations of initial state anisotropies, see discussion in Ref. [22|. 

These results are valid in the — 5 % centrality class as well. We have also checked that P(8v n ) is not sensitive 
to the freeze-out temperature within interval 100 < Tf Q < 160 MeV. Thus the distribution of relative fluctuations 
of v n could be the ideal observable to study the fluctuations of the initial geometry. If fluid dynamics provides a 
correct description of heavy-ion collisions, P(Sv n ) is a direct measurement of the initial state anisotropy fluctuations 
and can be compared directly to initial condition models, such as Monte-Carlo Glauber or various implementations 
of Color-Glass Condensate based initial conditions 1231. 



C. (vn,v m ) linear correlations 




FIG. 9. Correlations: a) c(v2,v 3 ), b) c(yz,Vi), and c) c{v2,Va), as function of transverse momentum in the 20 — 30 % centrality 
class using different initializations and viscosities. 
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c(v 3 ,v 4 ) 


sBC n/s = 0.0 


-0.09 


-0.11 


0.26 


0.32 


-0.03 


-0.11 


sBC rj/s = 0.16 


-0.09 


-0.11 


0.25 


0.63 


-0.03 


-0.09 


sWN rj/s = 0.16 


-0.15 


-0.14 


-0.04 


0.42 


0.03 


-0.11 



TABLE I. Linear correlation coefficients for p-r-integrated v n 's and e„'s in the 20 — 30 % centrality class. 

We computed the linear correlation coefficients, c(«2, ^3), c (t>2, Ui), and c(«3,W4) as a function of the transverse 
momentum, px- We found that c (v%, V3) ~ c (1)3, 1)4) ~ and, therefore, are not linearly correlated. These correlations 
are shown in Fig. and[3jD. We also show the values of the correlation coefficient between the integrated v n 's and 
of the coefficients e« in Table H] It should be noted that, even though both Glauber initializations used in this paper 
have the same relative anisotropy fluctuations, their correlations differ. 

Figure[9j; shows the correlation coefficient between V2 and V4 as a function of px- As can be read off from the Figure, 
c{v2,va) depends strongly on rj/s and, consequently, is sensitive to the properties of the QGP. Also, this correlation 
function is strongly affected by the correlation between €2 and £4 existing in the initial state, i.e., initial conditions 
with different 0(62,64) lead to very different 0(^2,^4), see the dashed and dotted lines in Fig. [Ht- Overall, just like V2 
probes the fluid-dynamical response to an initial geometry (characterized by £2), 0(^2,^4) probes the fluid-dynamical 
response to correlations in the initial geometry (in this case characterized by c(e2, £4))- 

In our investigations, c{v2,Vi) is the only correlation that was sensitive to both the fluctuations of the initial 
condition and the transport properties of the fluid. This correlation can thus be used as a further constraint to the 
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fluid-dynamical models applied to heavy-ion collisions. 



IV. CONCLUSIONS 



In this paper, we studied the relation between v n and e„ in ultrarelativistic heavy-ion collisions using event- by- event 
fluid dynamics. Wc confirmed that the second and third Fourier coefficients have a strong linear correlation to the 
initial geometry of the collision. We showed that while the event-average Fourier coefficients, {v n ) , n = 2, 3, and 4, 
are sensitive to the details of the fluid-dynamical evolution, their relative fluctuations, Sv n = (v n — (v n ) cv ) / (v n ) , 
are determined solely by the fluctuations of the corresponding initial state anisotropy coefficients, with basically no 
sensitivity to the viscosity of the fluid. This makes the distribution of 5v n a direct probe of the initial condition of 
a heavy-ion collision, providing a direct and clean measurement of the distribution of the relative fluctuations of the 
initial anisotropy, i.e., 

P (5v n ) ~ P (5e n ) , n = 2,3,4. (11) 

Surprisingly, this relation was shown to be true even for the relative fluctuations of V4, even though V4 itself is not 
linearly correlated to 64. 

Furthermore we found that the linear correlations between the flow harmonics v n are not solely dominated by the 
initial conditions. Especially the correlation function c(v2,Vi) is sensitive to the viscosity of the fluid providing an 
additional constraint for the model when one tries to extract the viscosity coefficient from the data. 
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Appendix: Relative widths of correlated distributions 



A relation between variables e and v (we omit the indexes n here) can be written as 

v = ke + d. 



If we define 



where c(e, v) is the linear correlator (Eq. 



k = —c{e,v), 



(A.l) 



(A.2) 



e and d are not linearly correlated, c(e, d) = 0. It also follows that 
a 2 d = <r 2 v (l-c(e,v) 2 ), (A3) 



i.e. the stronger the linear ( autocorrelation between e and d, the narrower the distribution of d. 

For the distributions of the scaled variables, Sx = (x — (x))/(x), it holds that ag x = a x /(x), and that (6x) = 0. 
The variance of Sv can now be written as 

1 



°5v = T~\2 [ k2(7 c + 2£;er £ cr rf c(e, d) + cr 2 ,] 



(v) 2 cr, 2 



(d) 
k(e) 



(v) 2 a 2 



(A.4) 



As mentioned, if c(e, v) w 1, w 0, and the last term in Eq. (|A.4[) is negligible. If also (d) ~ 0, <rs v w as e - Roughly 
speaking the requirement (d) w means that when e — > 0, (v) ~ — a requirement our distributions with n = 2 
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and 3 fulfill as seen in Figs. [TJ [2] [4j and [5l On the other hand, the correlation is weak for n — 4 in the 20 — 30 % 
centrality class. In that case the d distribution is wide, but (d) ^ as well, and the two terms in Eq. (|A.4[) sum to 
approximately one, and the distributions of 5e and 5v are equal even in that case. 
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